Cluster Explosive Synchronization in Complex Networks 
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The emergence of explosive synchronization has been reported as an abrupt transition in complex 
networks of first-order Kuramoto oscillators. In this Letter, we demonstrate that the nodes in a 
second-order Kuramoto model, perform a cascade of transitions toward a synchronous macroscopic 
state, which is a novel phenomenon that we call cluster explosive synchronization. We provide a 
rigorous analytical treatment using a mean-field analysis in uncorrelated networks. Our findings 
are in good agreement with numerical simulations and fundamentally deepen the understanding of 
microscopic mechanisms toward synchronization. 
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In the past few years, much research effort has been 
devoted to investigating the influence of network organi- 
zation on dynamical processes, such as random walks [1], 
congestion [H, Q , epidemic spreading ^] and synchroniza- 
tion d, m . Regarding synchronization of coupled oscil- 
lators, it has been demonstrated that the emergence of 
collective behavior in these structures depends on the 
patterns of connectivity of the underlying network. For 
instance, through a mean- field analysis, it has been found 
that Kuramoto oscillators display a second-order phase 
transition to the synchronous state with a critical cou- 
pling strength that depends on the network topology (Hj . 

Recently, discontinuous transitions to phase synchro- 
nization have been observed in SF networks |6|-|9|- 
This phenomenon, called explosive synchronization, was 
proved to be caused exclusively by a microscopic corre- 
lation between the network topology and the intrinsic 
dynamics of each oscillator. More specifically, Gomez- 
Gardefiez et al. Q considered the natural frequencies 
positively correlated with the degree distribution of the 
network, defining the natural frequency of each oscillator 
as equal to its number of connections. 

In this Letter, we substantially extend a first-order Ku- 
ramoto model used in Q to a second-order Kuramoto 
model fiol-fist that we modify in order to analyze global 
synchronization, considering the natural frequency of 
each node proportional to its degree In this 

model, we find a discontinuous phase transition in which 
small degree nodes join the synchronous component si- 
multaneously, whereas other nodes synchronize succes- 
sively according to their degrees (in contrast to @ where 
all the nodes join the synchronous component abruptly): 
This is a novel phenomenon which we call cluster explo- 
sive synchronization. By developing a mean-field theory 
we derive self-consistent equations that produce the lower 



and upper critical coupling strength associated to a hys- 
teretic behavior of the synchronization for uncorrelated 
networks. The analytical results are in good agreement 
with numerical simulations. Moreover, we show that de- 
creasing the network average frequency and increasing 
the coupling strength are the key factors that lead to 
cluster explosive synchronization. 

The second-order Kuramoto model consists of the fol- 
lowing set of equations fiol-list: 
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where Oi is the phase of the oscillator i = 1, . . . , A/", a 
the dissipation parameter, fli the the natural frequency 
distributed with a given probability density g{^), Xij the 
coupling strength and Aij an element of the network's 
adjacency matrix A. Here we assume that all connections 
have the same coupling which leads to a homogeneous 
coupling strength Xij = A, V z, j. To study the influence 
of dynamics and structure on global synchronization, we 
assume 
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where ki is the degree of node z, (k) the network aver- 
age degree and D a proportionality constant. A mean- 
field analysis allows us to investigate the dynamics of the 
model. 

We follow the continuum limit approach proposed 
in 01 and assume zero degree correlation between the 
nodes in the network. Denote by p{k] 0^ t) the frac- 
tion of nodes of degree k that have phase 6 at time 
t. The distribution p{k; 6>, t) is normalized according to 
J^^ p{k;6^t)d6 = 1. In an uncorrelated network, the 
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probability that a randomly picked edge is connected 
to a node with degree /c, phase 6 at time t is given 
by kP{k)p{k;0^t)/ (k)^ where P{k) is the degree dis- 
tribution and {k) the average degree. For a network, 
the order parameter r is defined as the magnitude of 
j,^iip{t) _ kiC^^^^'^^ I ki. In the continuum limit this 
yields 

re^^= / dk \ deP{k)k/{k)p{k',e,t)e'^^''\ (3) 

where kmin is the network's minimum degree and ij) is 
the average phase. When the phases are randomly dis- 
tributed, r ^ 0. On the other hand, when the oscillators 
evolve with similar phases, r ^ 1. 

Seeking to write a continuum-limit version of Eq. [T] 
with the natural frequency D{k—{k)) and the constant A 
in terms of the mean-field quantities r and we multiply 
both sides of Eq. [3] by e~'^^ and take the imaginary part, 
obtaining 

9 = -aO + D{k - (A:)) + kXr sin(V^ - 6>), (4) 
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FIG. 1: (Color online) Parameter space of the pendulum (Eq. 
[8]) . Plotted is the external constant torque / versus the damp- 
ing strength /3. The red area indicates parameter combina- 
tions that give rise to a globally stable fixed point. In the 
white area, only a stable limit cycle exists. Yellow indicates 
the area of bi- stability: both fixed point and limit cycle are 
stable. 



which is the same equation that describes the movement 
of a damped driven pendulum. 

In order to derive sufficient conditions for synchroniza- 
tion, we choose a reference frame that rotates with the 
average phase il) of the system, defining 0(t) = 0{t) — ip{t). 
Substituting the transformed variables (j){t) into the equa- 
tions of motion Eq. |4]and defining C{Xr) = {i/j -\- aip) / D , 
we get 

4) = -a(j) + D{k - (k) - C{Xr)) - kXr sin 0. (5) 

Eq. [5]can be interpreted as an extension to the second- 
order case of the model recently proposed in jij. The 
first-order Kuramoto model studied in jij presents hys- 
teresis only when SF topologies are considered in which 
each node's natural frequency is proportional to its de- 
gree. In contrast, it is known that systems described by 
a second-order Kuramoto model present hysteresis inde- 
pendently of the choice of the natural frequency distri- 
bution tii,[iii. 

In order to obtain sufficient conditions for the existence 
of the synchronous solution of Eq. [3l we derive a self- 
consistent equation for the order parameter r that can 
be written as the sum of the contribution riock due to 
oscillators which are phase-locked to the mean-field and 
the contribution of the non- locked drift oscillators r drift, 
i.e., r = riock ^ rdrift 0- 

Locked oscillators are characterized by ^ = ^ = 0, and 
turn out to possess degrees in a certain range k G [/ci, /C2]. 
Each locked oscillator has a /c-dependent constant phase. 



oscillators we obtain 
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independent single-peaked distribution. For the locked 
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On the other hand, the phase of the drift oscilla- 
tors rotates with period T in the stationary state, so 
that their density p{k](j)^t) satisfies p ~ \<P\~'^ [18]- As 
§ p{k](j))d(j) = Jq p{k](j))(j)dt = 1, this implies p{k;(j)) = 
T-^\(j)-^\ After substituting p{k](t)) into Eq.[3]and 
performing some mathematical manipulations motivated 
by p^], we get 
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The order parameter r is determined by the sum of Eq. [6] 
and [71 

It is known that systems subject to the equations of 
motion given by Eq. [H present a hysteresis as A is var- 
ied fisl , [2Q| . Therefore we consider in the following the 
system's dynamics for two distinct cases: (i) Increasing 
the coupling strength A. In this case, the system starts 
without synchrony (r ^ 0) and, as A is increased, ap- 
proaches the synchronous state (r ^ 1). (ii) Decreasing 
the coupling strength A. Now the system starts at the 
synchronous state (r ~ 1) and, as A is decreased, ever 
more oscillators lose synchronism, falling into the drift 
state. 

For the case in which the coupling strength A is in- 
creased from Ao, the synchronous state emerges after 
a threshold Xi has been crossed. Here we derive self- 
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consistent equations that allow to compute A^. In order 
to do that, we first investigate how the range of degree 
k of oscillators that are entrained in the mean-field de- 
pends on A. For convenience, we change the time-scale 
in Eq. [5]to r = V^Art, which yields 
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where we define P = a/VXkr and / = D{k — (k) — 
C{Xr))/kXr. The variable P is the damping strength 
and / corresponds to a constant torque (cf. the damped 
driven pendulum (2QI . [2i|). Our change of time-scale 
allows us to employ Melnikov's analysis [21] to deter- 
mine the range of integration [/ci,/^!] in the calculation 
of 
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Using Melnikov's analysis [18|, |20|, |2l|| we find that, for 
/ > 1, Eq. [8] has only one stable limit cycle solution 
(see Fig. [1]). If4/3/7r</<l, the system is bistable 
and a synchronized state co-exists with the limit cycle 
solution. If the coupling strength is increased further by 
a small amount, the synchronized state can only exist for 
I < 4/3/7r, where Eq. [8] has a stable fixed point solution, 
even for damping values close to /3 :^ 1 [li, Hlj . By 
solving the inequalities sine/) < 1 and / < 4/3/7r, we get 
the following range of k^ for the phase-locked oscillators 



k ^ [^1 7 



-fi — 



B - a/ 52 _ 4D4 + C{\r)Y 



2Z)2 



B + ^JB'^ - 4D4 {{k) + C{\r)y 
2D2 



(9) 



where 



B = 2D^{{k)^C{\r))^ 



If we now substitute Eq. [9] into the self-consistent Eqs. [6] 
and [71 we obtain and A^. 

When the coupling strength A is decreased, the oscilla- 
tors start at the phase- locked synchronous state, reach- 
ing the asynchronous state after a threshold A^. In 
order to calculate the threshold, we again investigate 
the range of degree k of phase- locked oscillators. Im- 
posing the phase locked solution in Eq. O we obtain 
sin0 = \D{k-{k)-C{Xr))\ ^ ^ ^^^^ ^^^^^^ 

oscillators are the nodes with degree k in the following 
range as a function of Ar 
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This allows us to calculate and A^ from the self- 
consistent Eqs. [6] and 

To check the validity of our mean- field analysis, we 
conduct numerical simulations of the model with a = 
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FIG. 2: (Color online) Comparison between numerical and 
analytical results. Shown is the order parameter for increasing 
(decreasing) A based on: i) simulations, red (blue) line and 
ii) self-consistent equations ( [6|) with synchronized degree 
/c^in Eq.IIOlA:^ in Eq. E)) , green (black) line. 



0.1 and D = 0.1 on SF networks characterized by TV = 
3000, (/c) = 10, kmin = 5 and the degree distribution 
P(/c) ~ k~^ with 7 = 3. Again, because we anticipate 
hysteresis, we have to distinguish two cases: first, we 
gradually increase A from Aq by amounts of (5 A, and for 
A = Ao, Ao + (5A, Ao + n(5A compute the increasing order 
parameter r^, where ^A = 0.1. Second, we gradually 
decrease A from Aq + nJA back toward Aq by amounts of 
(5A, this time computing the decreasing order parameter 
r^. Before each (5A-step, we simulate the system long 
enough to arrive at an attr actor. 

According to Fig. O the dependence of the simulated 
order parameter r on A indeed shows the expected hys- 
teresis: for increasing A, the discontinuous transition to 
synchrony happens at a A^ that is far larger than the 
threshold A^ of the backwards transition for decreasing 
A. Seeking to compare these simulation results to our 
mean- field analysis, we simultaneously solve Eqs. (0 El 
[9]) (resp. Eqs. ([HI [71 [TQ|) ) for the increasing (resp. de- 
creasing) case numerically to obtain analytical curves of 
r. 

A key ingredient of the solution process is C(Ar) which 
we retrieve from the simulation data. More specifically, 
recalling that C(Ar) just depends on and we assume 
that (i) C(Ar) ^ before the transition to synchrony, 
as the nodes oscillate independently and produce an un- 
changing zero mean field, and that (ii) C(Ar) ~ aip / D 
after the transition to synchrony, as the synchronized 
nodes dominate and produce a mean field that rotates 
with a constant frequency ij) (which we take from the 
simulations, cf. Fig. [3]). Fig. [2] reveals that our mean- 
field analysis predicts the critical thresholds A^ and A^ 
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FIG. 3: (Color online) Results from the numerical simulations 
with increasing coupling strength. The red line shows the 
mean value of '0(A), the red shade its standard deviation. The 
inset contains a periodic time series of at coupling strength 
0.8. 



very well. 

What happens at the node level when the transition 
to synchrony occurs? As just stated, the average fre- 
quency 2p that enters the equations via C{Xr) is a cru- 
cial quantity. Seeking to understand how the differ- 
ent nodes contribute to ip, for increasing A we calculate 
the average frequency of all nodes of degree /c, (cj)^ = 

E[.\k,=k]^^/iNPik)), where = J^+^ Ur)dt/T. As 
Fig. \^a) reveals, nodes of the same degree form clusters 
that join the synchronous component successively, start- 
ing from small degrees. This is in sharp contrast to the 
discontinuous phase transition observed in [6] , where the 
average frequency (cj)^ jumps to (k) at A^ for all k at 
the same time. We call this newly observed phenomenon 
cluster explosive synchronization. 

In Fig. Sl^b), we show the evolution of the lower and 
upper limits of the range of synchronous degrees, /c( and 
/c2, as a function of the coupling strength A. Analytical 
and simulation results are again in good agreement. Note 
the discontinuity in the evolution of /c( and /c^ which 
gives rise to the discontinuous transition in Fig. [2l For 
A > A^, the lower limit /c( is kept constant, and the higher 
limit of the synchronous nodes grows linearly with A 
(Fig.Hb)). 

In summary, we have demonstrated that a cluster ex- 
plosive synchronization transition occurs in a second- 
order Kuramoto model. As in previous studies on ex- 
plosive synchronization, a correlation between dynamics 
(natural frequency of a node) and local topology (the 
node's degree) is a necessary condition. With the con- 
nection between natural frequencies and local topology, 
we have presented the first analytical treatment of clus- 
ter explosive synchronization which is based on a mean- 
field approach in uncorrelated complex networks. Our 
simulations are in good agreement with the theory. Fur- 
thermore, we have shown that clusters of nodes of the 
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FIG. 4: (Color online) (a) Results from numerical simulations 
with increasing coupling strength A. Shown is the average fre- 
quency of nodes of various degrees against A. (b) synchronized 
degrees from simulations and mean field analysis. 



same degree join the synchronous component succes- 
sively, starting with small degrees. Our findings enhance 
the understanding of cluster explosive synchronization in 
the macrostate of a system and its applications will have 
a strong impact on the detection of clusters in larger 
networks. Also, our first analytical treatment can be ex- 
tended to applications |1QI413] where the use of second- 
order Kuramoto oscillators is relevant. 
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